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A previously-developed hybrid particle-continuum method [J. B. Bell, A. Garcia and S. 
A. Williams, SIAM Multiscale Modeling and Simulation, 6:1256-1280, 2008] is generalized 
to dense fluids and two and three dimensional flows. The scheme couples an explicit fluc- 
tuating compressible Navier-Stokes solver with the Isotropic Direct Simulation Monte Carlo 
(DSMC) particle method [A. Donev and A. L. Garcia and B. J. Alder, ArXiv preprint 
0908.0510]. To achieve bidirectional dynamic coupling between the particle (microscale) 
and continuum (macroscale) regions, the continuum solver provides state-based boundary 
conditions to the particle subdomain, while the particle solver provides flux-based bound- 
ary conditions for the continuum subdomain. This type of coupling ensures both state and 
flux continuity across the particle-continuum interface analogous to coupling approaches for 
deterministic parabolic partial differential equations; here, when fluctuations are included, 
a small (< 1%) mismatch is expected and observed in the mean density and temperature 
across the interface. By calculating the dynamic structure factor for both a "bulk" (periodic) 
and a finite system, it is verified that the hybrid algorithm accurately captures the prop- 
agation of spontaneous thermal fluctuations across the particle-continuum interface. The 
equilibrium diffusive (Brownian) motion of a large spherical bead suspended in a particle 
fluid is examined, demonstrating that the hybrid method correctly reproduces the velocity 
autocorrelation function of the bead but only if thermal fluctuations are included in the 
continuum solver. Finally, the hybrid is applied to the well-known adiabatic piston problem 
and it is found that the hybrid correctly reproduces the slow non-equilibrium relaxation of 
the piston toward thermodynamic equilibrium but, again, only the continuum solver includes 
stochastic (white-noise) flux terms. These examples clearly demonstrate the need to include 
fluctuations in continuum solvers employed in hybrid multiscale methods. 
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I. INTRODUCTION 

With the increased interest in nano- and micro-fluidics, it has become necessary to develop tools 
for hydrodynamic calculations at the atomistic scale [H El E]- While the Navier-Stokes-Fourier 
continuum equations have been surprisingly successful in modeling microscopic flows [3] , there are 
several issues present in microscopic flows that are difficult to account for in models relying on a 
purely PDE approximation. For example, it is well known that the Navier-Stokes equations fail 
to describe flows in the kinetic regions (large Knudsen number flows) that appear in small-scale 
gas flows [5j. It is also not a priori obvious how to account for the bidirectional coupling between 
the flow and embedded micro-geometry or complex boundaries. Furthermore, it is not trivial to 
include thermal fluctuations in Navier-Stokes solvers [6] [7] |9], and in fact, most of the time 
the fluctuations are omitted even though they can be important at instabilities [10] or in driving 
polymer dynamics |11[ [12]. An alternative is to use particle-based methods, which are explicit 
and unconditionally stable, robust, and simple to implement. The fluid particles can be directly 
coupled to the microgeometry, for example, they can directly interact with the beads of a polymer 
chain. Fluctuations are naturally present and can be tuned to have the correct spatio-temporal 
correlations. 

Several particle methods have been described in the literature, such as molecular dynamics 
(MD) P3], Direct Simulation Monte Carlo (DSMC) p3], dissipative particle dynamics (DPD) [15] . 
and multi-particle collision dynamics (MPCD) |16|ll7j. Here we use the Isotropic DSMC (I-DSMC) 
stochastic particle method described in Ref. [18] . In the I-DSMC method, deterministic interactions 
between the particles are replaced with stochastic momentum exchange (collisions) between nearby 
particles. The I-DSMC method preserves the essential hydrodynamic properties of expensive MD: 
local momentum conservation, linear momentum exchange on length scales comparable to the 
particle size, and a similar fluctuation spectrum. At the same time, the I-DSMC fluid is ideal and 
structureless, and as such is much simpler to couple to a continuum solver. 

However, even particle methods with coarse-grained dynamics, such as I-DSMC, lack the effi- 
ciency necessary to study realistic problems because of the very large numbers of particles needed to 
fill the required computational domain. Most of the computational effort in such a particle method 
would, however, be expended on particles far away from the region of interest, where a description 
based on the Navier-Stokes equations is likely to be adequate. Hybrid methods are a natural can- 
didate to combine the best features of the particle and continuum descriptions. A particle method 
can be used in regions where the continuum description fails or is difficult to implement, such as 
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near suspended structures or complex boundaries, and a more efficient continuum description can 
be used to around the particle domain, as illustrated in Fig. [T] This type of hybrid algorithm fits 
in the Multi-Algorithm Refinement (MAR) simulation approach to the modeling and simulation 
of multiscale problems |1 9|, 120) [21]. MAR combines two or more algorithms, each of which is ap- 
propriate for a different scale regime. The general idea is to perform detailed calculations using 
an accurate but expensive (e.g., particle) algorithm in a small region, and couple this computa- 
tion to a simpler, less expensive (e.g., continuum) method applied to the rest. The challenge is 
to ensure that the numerical coupling between the different levels of the algorithm hierarchy, and 
especially the coupling of the particle and continuum computations, is self-consistent, stable, and 
most importantly, does not adversely impact the underlying physics. 




Figure 1 : A two-dimensional illustration of the use of the MAR hybrid to study a polymer chain (larger red 
circles) suspended in an I-DSMC particle fluid. The region around the chain is filled with particles (smaller 
green circles), while the remainder is handled using a fluctuating hydrodynamic solver. The continuum 
(macro) solver grid is shown (thicker blue lines), along with the (micro) grid used by the particle method 
(thinner blue lines) . The fluctuating velocities in the continuum region are shown as vectors originating from 
the center of the corresponding cell (purple). The interface between the particle and continuum regions is 
highlighted (thicker red line). 

A crucial feature of our hybrid algorithm is that the continuum solver includes thermal fluctu- 
ations in the hydrodynamic equations consistent with the particle dynamics, as previously investi- 
gated in one dimension in Ref. [2T] . Thermal fluctuations play an important role in describing the 
state of the fluid at microscopic and mesoscopic scales, especially when investigating systems where 
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the microscopic fluctuations drive a macroscopic phenomenon such as the evolution of instabilities, 
or where the thermal fluctuations drive the motion of suspended microscopic objects in complex 
fluids. Some examples in which spontaneous fluctuations can significantly affect the dynamics in- 
clude the breakup of droplets in jets |22[ [23] , Brownian molecular motors [H] , Rayleigh-Bernard 
convection [25] , Kolmogorov flows |26[ [27] , Rayleigh- Taylor mixing [10] , and combustion and ex- 
plosive detonation [2H] • In our algorithm, the continuum solver is a recently-developed three-stage 
Runge-Kutta integration scheme for the Landau-Lifshitz Navier-Stokes (LLNS) equations of fluc- 
tuating hydrodynamics in three dimensions [9], although other finite- volume explicit schemes can 
trivially be substituted. 

As summarized in Section [TTJ the proposed hybrid algorithm is based on a fully dynamic bidirec- 
tional state-flux coupling between the particle and continuum regions. In this coupling scheme the 
continuum method provides state-based boundary conditions to the particle subdomain through 
reservoir particles inserted at the boundary of the particle region at every particle time step. During 
each continuum time step a certain number of particle time steps are taken and the total particle 
flux through the particle-continuum interface is recorded. This flux is then imposed as a flux-based 



boundary condition on the continuum solver, ensuring strict conservation |20[ 12 lj . Section III 



describes the technical details of the hybrid algorithm, focusing on components that are distinct 
from those described in Refs. |20[ 12 lj . Notably, the use of the Isotropic DSMC particle method 
instead of the traditional DSMC method requires accounting for the interactions among particles 
that are in different continuum cells. 



In Section IV we thoroughly test the hybrid scheme in both equilibrium and non-equilibrium 
situations, and in both two and three dimensions. In Section IV A| we study the continuity of 
density and temperature across the particle-continuum interface and identify a small mismatch of 
order 1/Nq, where iVo is the number of particles per continuum cell, that can be attributed to the 



use of fluctuating values instead of means. In Section IV B we compute dynamic structure factors 
in periodic ("bulk") and finite quasi two- and one-dimensional systems and find that the hybrid 
method seamlessly propagates thermal fluctuations across the particle-continuum interface. 

In Section IV C| we study the diffusive motion of a large spherical buoyant bead suspended in 
a bead of I-DSMC particles in three dimensions by placing a mobile particle region around the 
suspended bead. This example is of fundamental importance in complex fluids and micro-fluidics, 
where the motion of suspended objects such as colloidal particles or polymer chains has to be 
simulated. Fluctuations play a critical role since they are responsible for the diffusive motion of the 
bead. The velocity-autocorrelation function (VACF) of a diffusing bead has a well-known power 
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law tail of hydro dynamic origin and its integral determines the diffusion coefficient. Therefore, 
computing the VACF is an excellent test for the ability of the hybrid method to capture the 
influence of hydrodynamics on the macroscopic properties of complex fluids. 



Finally, in Section IV D we study the slow relaxation toward thermal equilibrium of an adiabatic 
piston with the particle region localized around the piston. In the formulation that we consider, the 
system is bounded by adiabatic walls on each end and is divided into two chambers by a mobile and 
thermally insulating piston that can move without friction. We focus on the case when the initial 
state of the system is in mechanical equilibrium but not thermodynamic equilibrium: the pressures 
on the two sides are equal but the temperatures are not. Here we study the slow equilibration of 
the piston towards the state of thermodynamic equilibrium, which happens because asymmetric 
fluctuations on the two sides of the piston slowly transfer energy from the hotter to the colder 
chamber. 

We access the performance of the hybrid by comparing to purely particle simulations, which 
are assumed to be "correct". Unlike in particle methods, in continuum methods we can trivially 
turn off fluctuations by not including stochastic fluxes in the Navier-Stokes equations. By turning 
off fluctuations in the continuum region we obtain a deterministic hybrid scheme, to be contrasted 
with the stochastic hybrid scheme in which fluctuations in the continuum region are consistent with 
those in the particle region. By comparing results between the deterministic and stochastic hybrid 
we are able to assess the importance of fluctuations. We find that the deterministic hybrid gives the 
wrong long-time behavior for both the diffusing spherical bead and the adiabatic piston, while the 
stochastic hybrid correctly reproduces the purely particle runs at a fraction of the computational 
cost. These examples demonstrate the need to include thermal fluctuations in the continuum 
solvers in hybrid particle-continuum methods. 



II. BRIEF OVERVIEW OF THE HYBRID METHOD 



In this section we briefly introduce the basic concepts behind the hybrid method, delegating 
further technical details to later sections. Our scheme is based on an Adaptive Mesh and Algo- 
rithm Refinement (AMAR) methodology developed over the last decade in a series of works in 
which a DSMC particle fluid was first coupled to a deterministic compressible Navier-Stokes solver 
in three dimensions |20| [29] , and then to a stochastic (fluctuating) continuum solver in one dimen- 
sion [21]. This section presents the additional modifications to the previous algorithms necessary 
to replace the traditional DSMC particle method with the isotropic DSMC method [18], and a 
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full three-dimensional dynamic coupling of a complex particle fluid [30] with a robust fluctuating 



hydrodynamic solver [9]. These novel techniques are discussed further in Section III 



Next, we briefly describe the two components of the hybrid, namely, the particle microscopic 
model and the continuum macroscopic solver, and then outline the domain decomposition used to 
couple the two, including a comparison with other proposed schemes. Both the particle algorithm 
and the macroscopic solver have already been described in detail in the literature, and furthermore, 
both can easily be replaced by other methods. Specifically, any variant of DSMC and MPCD can 
be used as a particle algorithm, and any explicit finite-volume method can be used as a continuum 
solver. For this reason, in this paper we focus on the coupling algorithm. 



A. Particle Model 



The particle method that we employ is the Isotropic DSMC (I-DSMC) method, a dense fluid 1 
generalization of the Direct Simulation Monte Carlo (DSMC) algorithm for rarefied gas flows p3] . 
The I-DSMC method is described in detail in Ref. [18] and here we only briefly summarize it. It is 
important to note that, like the traditional DSMC fluid, the I-DSMC fluid is an ideal fluid, that is, 
it has the equation of state (EOS) and structure of an ideal gas; it can be thought of as a viscous 
ideal gas. As we will see shortly, the lack of structure in the I-DSMC fluid significantly simplifies 
coupling to a continuum solver while retaining many of the salient features of a dense fluid. 

In the I-DSMC method, the effect of interatomic forces is replaced by stochastic collisions 
between nearby particles. The interaction range is controlled via the collision diameter D or, 
equivalently, the density (hard-sphere volume fraction) (ft = irND^/^QV), where N is the total 
number of particles in the simulation volume V. The strength of the interaction is controlled 
through the dimensionless cross-section pref actor x, which is on the order of unity [18] . Stochastic 
collisions are processed at the end of every particle time step of duration Atp, and in-between 
collisions each particle i is streamed advectively with constant velocity Vi = f"j. The spatial 
domain of the simulation, typically a rectangular region, is divided into micro cells of length 
L c ft D, which are used to efficiently identify all particles that are within distance L c of a given 
particle by searching among the particles in neighboring cells (each cell has 3 d neighboring cells, 
counting itself, where d is the spatial dimension). At each time step, conservative collisions occur 
between randomly-chosen pairs of particles that are closer than a distance D apart; specifically, 



1 Note that by a dense fluid we mean a fluid where the mean free path is small compared to the typical fluid 
inter-atomic distance. 
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a random amount of momentum and energy is exchanged between the two particles with the 
probabilities of the collision outcomes obeying detailed balance. The I-DSMC algorithm can be 
viewed as a stochastic alternative to deterministic hard-sphere molecular dynamics (MD) [31 , 
where hard spheres of diameter D collide when they touch. 

In addition to the fluid (I-DSMC) particles there may be a number of additional spherical 
particles, which we refer to as beads, suspended in the I-DSMC fluid. The beads interact with 
each other and the fluid particles either deterministically, as hard spheres impermeable to and 
colliding with touching particles, or stochastically, as permeable spheres colliding stochastically with 
overlapping particles. Note that the sphere radii used for determining the fluid-fluid, fluid-bead, 
and bead-bead interaction distances may, in general, be different. Combining deterministic hard- 
sphere collisions with stochastic collisions requires a mixed time-driven and event-driven treatment, 
as in the Stochastic Event-Driven MD (SEDMD) algorithm developed in Ref. (30]. We will use 
the SEDMD algorithm in the examples presented in this paper. 

We have also developed an alternative purely time-driven fluid-bead coupling in which the fluid 
is allowed to permeate the beads and all of the particle interactions are stochastic. This leads to 
a much simpler algorithm that can also easily be parallelized. This distinction between hard and 
permeable spheres has analogues in other methods for complex fluids. For example, in Lattice- 
Boltzmann simulations (32], beads can either be modeled as hard spheres (using a bounce-back 
collision rule at the lattice sites on the surface of the bead), or, more efficiently and commonly, as 
permeable spheres that let the fluid pass through them but experience a frictional force due to the 
fluid motion (exerting the opposite force back on the lattice sites they overlap with). 



B. Continuum Model 



At length scales and time scales larger than the molecular ones, the dynamics of the particle 
fluid can be coarse grained |33|[3"3] to obtain evolution equations for the slow macroscopic variables. 
Specifically, we consider the continuum conserved fields 



U(r,t) 



U(r,t) = J2 
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(1) 



where the conserved variables, namely the densities of mass p, momentum j = pv, and energy 
e = e(p, T) + \pv 2 , can be expressed in terms of the primitive variables, mass density p, velocity v 
and temperature T; here e is the internal energy density. Here the symbol = means that we consider 
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a stochastic field U(r, t) that approximates the behavior of the true atomistic configuration U(r, t) 
over long length and time scales (compared to atomistic scales) in a certain integral average sense; 
notably, for sufficiently large cells the integral of U(r,t) over the cell corresponds to the total 
particle mass, momentum and kinetic energy contained inside the cell. 

The evolution of the field U(r, t) is modeled with the Landau-Lifshitz Navier-Stokes (LLNS) 
system of stochastic partial differential equations (SPDEs) in d dimensions, given in conservative 
form by 



d t U = -V ■ [F(U) - Z(U,r,t)] 



(2) 



where the deterministic flux is taken from the traditional compressible Navier-Stokes-Fourier equa- 
tions, 



F(U) 



pv 



pvv T + PI — a 
(e + P)v - {a ■ v + £) 

where P = P(p, T) is the pressure, the viscous stress tensor is a = 2r] 
(we have assumed zero bulk viscosity), and the heat flux is £ = /iVT. As postulated by Landau- 
Lifshitz [Ml ES] , the stochastic flux 

r 
S 

£ • v + H 

is composed of the stochastic stress tensor I] and stochastic heat flux vector H, assumed to be 
mutually uncorrelated random Gaussian fields with a covariance 

2 



Z = 



<E(r, t)E*(r' s t')) =C^5{t - t')5{r - r'), where cg£ = 2fjk B T ( fofy + 5 u S jk - 



[S{r,t)S*(r',t')) =C s 5{t - t')8(r - r'), where cf S) 



2pk B T5ij, 



(3) 



where overbars denote mean values. 

As discussed in Ref. [9], the LLNS equations do not quite make sense written as a system of 
nonlinear SPDEs, however, they can be linearized to obtain a well-defined linear system whose 
equilibrium solutions are Gaussian fields with known covariances. We use a finite-volume dis- 
cretization, in which space is discretized into N c identical macro cells Vj of volume V c , and the 
value Uj stored in cell 1 < j < N c is the average of the corresponding variable over the cell 



Uj(t) 



1 



U(r,t)dr 
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c jv 



U(r,t)dr, 
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where U is defined in Eq. Time is discretized with a time step Ate, approximating U(r,t) 
pointwise in time with U n = {J7", U%^\, 

where n > enumerates the macroscopic time steps. While not strictly necessary, we will assume 
that each macro cell consists of an integer number of micro cells (along each dimension of the 
grid), and similarly each macro time step consists of an integer number n ex of micro time steps, 
At c = n ex At P . 

In addition to the cell averages Uj", the continuum solver needs to store the continuum normal 
flux Fj'j, through each interface I = Vj Pi Vy between touching macro cells j and f during a given 
time step, 

3 

where Sjji is the surface area of the interface, and F^ji = ~Ff,j- Here we will absorb the 
various prefactors into a total transport (surface and time integrated flux) through a given macro 
cell interface /, = V~ 1 AtSjF r }, which simply measures the total mass, momentum and energy 
transported through the surface I during the time interval from time t to time t+At. We arbitrarily 
assign one of the two possible orientations (direction of the normal vector) for each cell-cell interface. 
How the (integrated) fluxes 3?" are calculated from U" does not formally matter; all that the hybrid 
method uses to advance the solution for one macro time step are C7™, U™ +1 and . Therefore, any 
explicit conservative finite-volume method can be substituted trivially. Given this generality, we 
do not describe in any detail the numerical method used to integrate the LLNS equations; readers 
can consult Ref. [9] for further information. 



C. Coupling between particle and continuum subdomains 

The hybrid method we use is based on domain decomposition, and is inspired by Adaptive Mesh 
Refinement (AMR) methodology for conservation laws |36 | l37j. Our coupling scheme closely follows 
previously-developed methodology for coupling a traditional DSMC gas to a continuum fluid, first 
proposed in the deterministic setting in Ref. [20] and then extended to a fluctuating continuum 
method in Ref. |21j . The key new ingredient is the special handling of the collisional momentum 
and energy transport across cell interfaces, not found in traditional DSMC. For completeness, we 
describe the coupling algorithm in detail, including components already described in the literature. 
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We split the whole computational domain into particle and continuum subdomains, which com- 
municate with each other through information near the particle-continuum interface I, assumed 
here to be oriented such that the flux measures the transport of conserved quantities from 
the particle to the continuum regions. In AMR implementations subdomains are usually logically 
rectangular patches; in our implementation, we simply label each macro cell as either a particle 
cell or a continuum cell based on whatever criterion is appropriate, without any further restrictions 
on the shape or number of the resulting subdomains. For complex fluids applications, macro cells 
near beads (suspended solute) and sometimes near complex boundaries will be labeled as particle 
cells. The continuum solver is completely oblivious to what happens inside the particle subdomain 
and thus it need not know how to deal with complex moving boundaries and suspended objects. 
Instead, the continuum solution feels the influence of boundaries and beads through its coupling 
with the particle subdomains. 

The dynamic coupling between particle and continuum subdomains is best viewed as a mu- 
tual exchange of boundary conditions between the two regions. Broadly speaking, domain- 
decomposition coupling schemes can be categorized based on the type of boundary conditions 
each subdomain specifies for the other [38]. Our scheme is closest to a state-flux coupling scheme 
based on the classification proposed in Ref. [38] for incompressible solvers (the term "velocity- 
flux" is used there since velocity is the only state variable). A state-flux coupling scheme is one 
in which the continuum solver provides to the particle solver the conserved variables U in the 
continuum reservoir macro cells near the particle-continuum interface /, that is, the continuum 
state is imposed as a boundary condition on the particle region. The particle solver provides to 
the continuum solver the flux through the interface /, that is, the particle flux is imposed as 
a boundary condition on the continuum subdomain. This aims to achieve continuity of both state 
variables and fluxes across the interface, and ensures strict conservation, thus making the coupling 
rather robust. Note that state/flux information is only exchanged between the continuum and 
particle subdomains every n ex particle (micro) time steps, at the beginning/end of a macro time 



step. A more detailed description of the algorithm is given in Section III 



1. Comparison with other coupling schemes 

There are several hybrid methods in the literature coupling a particle method, and in partic- 
ular, molecular dynamics (MD), with a continuum fluid solver [39 . There are two main types 
of applications of such hybrids [ID]. The first type are problems where the particle description 
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is localized to a region of space where the continuum description fails, such as, for example, a 
complex boundary, flow near a corner, a contact line, a drop pinchoff region, etc.. The second type 
are problems where the continuum method needs some transport coefficients, e.g., stress-strain 
relations, that are not known a priori and are obtained via localized MD computations. In the 
majority of existing methods a stationary solution is sought |41j . and a deterministic incompress- 
ible or isothermal formulation of the Navier-Stokes equations is used in the continuum |4U| |4*2] . By 
contrast, we are interested in a fully dynamic bidirectional coupling capable of capturing the full 
range of hydrodynamic effects including sound and energy transport. We also wish to minimize 
the size of the particle regions and only localize the particle computations near suspended objects, 
making it important to minimize the artifacts at the interface. As we will demonstrate in this 
work, including thermal fluctuations in the continuum formulation is necessary to obtain a proper 
coupling under these demanding constraints. 

The only other work we are aware of that develops a coupling between a fluctuating compressible 
continuum solver and a particle method, specifically molecular dynamics, is a coupling scheme 
developed over the last several years by Coveney, Flekkoy, de Fabritiis, Delgado-Buscallioni and 
collaborators [21 331 03]. There are two important differences between their method and our 
algorithm. Firstly, their method is (primarily, but not entirely) a flux-flux coupling scheme, unlike 
our state-flux coupling. Secondly, we do not use MD but rather I-DSMC, which, as discussed 
below, significantly simplifies the handling of the continuum-particle interface. 

It is not difficult to impose boundary conditions for the continuum subdomain based on the 
particle data, and any consistent boundary condition for the PDE being solved can in principle 
be imposed. When the continuum solver is deterministic the fluctuations (often inappropriately 
referred to as "noise") in the particle data need to be filtered using some sort of spatial and temporal 
averaging jUJ , or the continuum solver needs to be robust |20| 129] . 

The difficult part in coupling schemes is the handling of the boundary conditions for the particle 
subdomain. It is very difficult to truncate a fluid particle region without introducing artifacts in the 
structure of the particle fluid, and the transport coefficients of the particle fluid (e.g., the pressure 
and viscosity) critically depend on the fluid structure. In most molecular simulations periodic 
boundary conditions are used to avoid such artifacts, however, this is not possible in our case. 
In order to minimize artifacts at the boundary of the particle subdomain, most coupling schemes 
add an overlap region in addition to the reservoir region that we described above. In the overlap 
region, particles are simulated even though the region belongs to the continuum subdomain. The 
structure of the fluid in the overlap region is left to adjust to that in the particle subdomain, 
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thus minimizing the artifacts. At the same time, the dynamics of the particles in the overlap 
region is somehow constrained to match the underlying continuum subdomain dynamics using, for 
example, various forms of constrained Langevin-type thermostats or (non-holonomic) constrained 
MD [42, 45 . To prevent particles from flying outside of the particle subdomain, various artificial 
constraining forces are added in the reservoir region, and typically some particle insertion/deletion 
scheme is used as well. The details can be very different and tricky, and we are not aware of any 
detailed comparison or even rudimentary mathematical analysis of the different types of coupling 
other than the stability analysis of four idealized coupling schemes presented in Ref. [38J. 

We do not use an overlap region in the coupling scheme we present here. Firstly, and most 
importantly, because the I-DSMC fluid is structureless, an overlap region and the associated al- 
gorithmic complications are simply not required. The addition of an overlap region (added fluid 
mass) almost always introduces some delay and errors in the coupling, and it is usually assumed 
that this surface effect is negligible when the size of the overlap region is small compared to the par- 
ticle region and when the hydro time steps are much larger than the particle time step. These are 
assumptions that we do not wish to make as we try to minimize the size of the particle subdomain. 

The simple and direct coupling scheme we have presented, however, does not work if a structured 
stochastic particle fluid is used, such as, for example, a structured stochastic fluid like the Stochastic 
Hard-Sphere Dynamics (SHSD) fluid |18j . The SHSD fluid is weakly structured, so that at least 
particle insertion/deletion in the reservoir region is not a problem. However, an overlap region is 
necessary to smoothly match the fluid structure at the particle-continuum interface. Constraining 
the dynamics in the overlap region can be done in I-DSMC by introducing additional one-particle 
collisions (dissipation) that scatter the particle velocities so that their mean equals the continuum 
field values. However, it is not clear how to do this consistently when fluctuations are included 
in the continuum solver, and in this paper we restrict our focus on structureless (ideal) stochastic 
fluids. 



III. DETAILS OF THE COUPLING ALGORITHM 



The basic ideas behind the state-flux coupling were already described in Section II C In this sec- 
tion we describe in detail the two components of the particle-continuum coupling method, namely, 
the imposition of the continuum state as a boundary condition for the particle subdomain and the 
imposition of the particle flux as a boundary condition for the continuum subdomain. At the same 
time, we will make clear that our coupling is not purely of the state-flux form. The handling of 
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the continuum subdomain is essentially unchanged from the pure continuum case, with the only 
difference being the inclusion of a refluxing step. The handling of the particle subdomain is more 
complex and explained in greater detail, including pseudocodes for several steps involved in taking 
a micro time step, including insertion of reservoir particles and the tracking of the particle fluxes. 



A. State exchange 

We first explain how the state in the reservoir macro cells bordering the particle subdomain, 
denoted by U H , is used by the particle algorithm. The micro cells that are inside the reservoir 
macro cells and are sufficiently close to the particle subdomain to affect it during a time-interval 
Atp are labeled as reservoir micro cells. For I-DSMC fluids, assuming the length of the micro cells 
along each dimension is L c <J D, the micro cells immediately bordering the particle subdomain 
as well as all of their neighboring micro cells need to be included in the reservoir region. An 
illustration of the particle and reservoir regions is given in Figs. [T]and|2j 

At the beginning of each particle time step reservoir particles are inserted randomly into the 
reservoir micro cells. The number of particles inserted is based on the target density in the corre- 
sponding reservoir macro cell. The velocities of the particles are chosen from a Maxwell-Boltzmann 



or Chapman-Enskog [46 distribution (see discussion in Section III D 2 ) with mean velocity and tem- 
perature, and also their gradients if the Chapman-Enskog distribution is used, chosen to match the 
momentum and energy densities in the corresponding macro cell. The positions of the reservoir 
particles are chosen randomly uniformly (i.e., sampled from a Poisson spatial distribution) inside 
the reservoir cells which does not introduce any artifacts in the fluid structure next to the interface 
because the I-DSMC fluid is ideal and structureless. This is a major advantage of ideal fluids over 
more realistic structured non-ideal fluids. Note that in Ref. [30] we used a particle reservoir to 
implement open boundary conditions in pure particle simulations, the difference here is that the 
state in the reservoir cells comes from a continuum solver instead of a pre-specified stationary flow 
solution. 

The reservoir particles are treated just like the rest of the particles for the duration of a particle 
time step. First they are advectively propagated (streamed) along with all the other particles, and 
the total mass, momentum and energy transported by particles advectively through the particle- 
continuum interface, 4>p , is recorded. Particles that, at the end of the time step, are not in either 
a reservoir micro cell or in the particle subdomain are discarded, and then stochastic collisions are 
processed between the remaining particles. In the traditional DSMC algorithm, collisions occur 



14 




Figure 2: Illustration of a hybrid simulation of two-dimensional plug flow around a permeable stationary 
disk (red). The macroscopic grid is shown (dark-blue lines), with each macro cell composed of 6 x 6 micro 
cells (not shown) . The particle subdomain surrounds the disk and the particle-continuum interface is shown 
(red). A snapshot of the I-DSMC particles is also shown (green), including the reservoir particles outside of 
the particle subdomain. The time-averaged velocity in each continuum cell is shown, revealing the familiar 
plug flow velocity field that is smooth across the interface. This example clearly demonstrates that the 
continuum solver feels the stationary disk through the particle subdomain even though the continuum solver 
is completely oblivious to the existence of the disk. 

only between particles inside the same micro cell, and thus all of the particles outside the particle 
subdomain can be discarded [20, [21] . However, in the I-DSMC algorithm particles in neighboring 
cells may also collide, and thus the reservoir particles must be kept until the end of the time step. 
Collisions between a particle in the particle subdomain and a particle in a reservoir cell lead to 
collisional exchange of momentum and energy through the interface as well and this contribution 
must also be included in 3>p^. Note that the reservoir particles at the very edge of the reservoir 
region do not have an isotropic particle environment and thus there are artifacts in the collisions 
which they experience, however, this does not matter since it is only essential that the particles in 
the particle subdomain not feel the presence of the interface. 
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B. Flux exchange 

After the particle subdomain is advanced for n ex (micro) time steps, the particle flux <&p^ is 
imposed as a boundary condition in the continuum solver so that it can complete its (macro) 
time step. This flux exchange ensures strict conservation, which is essential for long-time stability. 
Assume that the continuum solver is a one-step explicit method that uses only stencils of width 
one, that is, calculating the flux for a given macro cell-cell interface only uses the values in the 
adjacent cells. Under such a scenario, the continuum solver only needs to calculate fluxes for the 
cell-cell interfaces between continuum cells, and once the particle flux <&p^ is known the continuum 
time step ^ can be completed. It is obvious that in this simplified scenario the coupling is purely 
of the state-flux form. However, in practice we use a method that combines pieces of state and flux 
exchange between the particle and continuum regions. 

First, the particle state is partly used to advance the continuum solver to the next time step. 
Our continuum solver is a multi-stage method and uses stencils of width two in each stage, thus 
using an effective stencils that can be significantly larger than one cell wide [9|. While one can 
imagine modifying the continuum solver to use specialized boundary stencils near the particle- 
continuum interface (e.g., one-sided differencing or extrapolation), this is not only more complex 
to implement but it is also less accurate. Instead, the continuum method solves for the fields over 
the whole computational domain (continuum patch in AMAR terminology) and uses hydrodynamic 
values for the particle subdomain (particle patch) obtained from the particle solver. These values 
are then used to calculate provisional fluxes <&r[ and take a provisional time step, as if there 
were no particle subdomain. This makes the implementation of the continuum solver essentially 
oblivious to the existence of the particle regions, however, it does require the particle solver to 
provide reasonable conserved values for all of the macro particle cells. This may not be possible 
for cells where the continuum hydrodynamic description itself breaks down, for example, cells that 
overlap with or are completely covered by impermeable beads or features of a complex boundary. 
In such empty cells the best that can be done is to provide reasonable hydrodynamic values, for 
example, values based on the steady state compatible with the specified macroscopic boundary 
conditions. For partially empty cells, that is, macro cells that are only partially obscured, one 
can use the uncovered fraction of the hydro cell to estimate hydrodynamic values for the whole 
cell. In practice, we have found that as long as empty cells are sufficiently far from the particle- 
continuum interface (in particular, empty cells must not border the continuum subdomain) the 
exact improvised hydrodynamic values do not matter much. Note that for permeable beads there 
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is no problem with empty cells since the fluid covers the whole domain. 

Secondly, the provisional continuum fluxes are partly used to advance the particle subdomain 
to the same point in time as the continuum solver. Specifically, a linear interpolation between the 
current continuum state and the provisional state is used as a boundary condition for the reservoir 
particles. This temporal interpolation is expected to improve the temporal accuracy of the coupling, 
although we are not aware of any detailed analysis. Note that if n ex = 1 this interpolation makes 
no difference and the provisional continuum fluxes are never actually used. 

Once the particle solver advances n ex time steps, a particle flux 3>p^ is available and it is 
imposed in the continuum solver to finalize the provisional time step. Specifically, hydrodynamic 
values in the particle macro cells are overwritten based on the actual particle state, ignoring the 
provisional prediction. In order to correct the provisional fluxes, a refluxing procedure is used in 
which the state U H in each of the continuum cells that border the particle-continuum interface 
are changed to reflect the particle 3>p^ rather than the provisional flux , 

jj{B) V (B) _ (I) (J) 

This refluxing step ensures strict conservation and ensures continuity of the fluxes across the 
interface, in addition to continuity of the state. 



C. Taking a macro time step 

Algorithm [T] summarizes the hybrid algorithm and the steps involved in advancing both the 
simulation time by one macro time step Ate = n ex Atp. Note that at the beginning of the simula- 
tion, we initialize the hydrodynamic values for the continuum solver, consistent with the particle 
data in the particle subdomain and generated randomly from the known (Gaussian) equilibrium 
distributions in the continuum subdomain. 

Algorithm 1: Take a macro time step by updating the continuum state U h from time t to time 

t + At c . 

1. Provisionally advance the continuum solver: Compute a provisional macro solution U^f xt at time 
t + Atc everywhere, including the particle subdomain, with an estimated (integrated) provisional 
flux <&£p through the particle-continuum interface. Reset the particle flux ffrp' 1 *— 0. 

2. Advance the particle solver: Take n ex micro time steps (see Algorithm [2]) : 

(a) At the beginning of each particle time step reservoir particles are inserted at the boundary of 
the particle subdomain with positions and velocities based on a linear interpolation between 
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U h and U^ xt (see Algorithm |4| . This is how the continuum state is imposed as a boundary 
condition on the particle subdomain. 

(b) All particles are propagated advectively by Atp and stochastic collisions are processed, accu- 
mulating a particle flux ffrp- 1 (see Algorithm ^ . 

3. Synchronize the continuum and particle solutions: 

(a) Advance: Accept the provisional macro state, U h <— U^ xt . 

(b) Correct: The continuum solution in the particle subdomain Up is replaced with cell averages 
of the particle state at time t + Ate, thus forming a composite state over the whole domain. 

(c) Reflux: The continuum solution U H in the macro cells bordering the particle subdomain is 
corrected based on the particle flux, <— — $>jp + $>p\ This effectively imposes the 
particle flux as a boundary condition on the continuum and ensures conservation in the hybrid 
update. 

(d) Update the partitioning between particle and continuum cells if necessary. Note that this step 
may convert a continuum cell into an unfilled (devoid of particles) particle cell. 



D. Taking a micro time step 

Taking a micro time step is described in Algorithm [2j and the remaining subsections give further 
details on the two most important procedures used. The initial particle configuration can most 
easily be generated by marking all macro cells in the particle domain as unfilled. 

Algorithm 2: Take an I-DSMC time step. We do not include details about the handling of the 
non-DSMC (solute) particles here. Note that the micro time step counter np should be 
re- initialized to zero after every n ex time steps. 

1. Visit all macro cells that overlap the reservoir region or are unfilled (i.e., recently converted from 
continuum to particle) one by one, and insert trial particles in each of them based on the continuum 
state U h, as described in Algorithm [4] 

2. Update the clock t <— t + Atp and advance the particle subdomain step counter np <— np + 1. Note 
that when an event-driven algorithm is used this may involve processing any number of events that 
occur over the time interval Atp [30j . 

3. Move all I-DSMC particles to the present time, updating the total kinetic mass, momentum and 
energy flux through the coupling interface F^p whenever a particle crosses from the particle into the 
continuum subdomain or vice versa, as detailed in Algorithm [3] 
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4. Perform stochastic collisions between fluid particles, including the particles in the reservoir region. 
Keep track of the total collisional momentum and energy flux through the coupling interface by 
accounting for the amount of momentum Ap^ = mAvij and kinetic energy Ae^ transferred from a 
regular particle i (in the particle subdomain) to a reservoir particle j (in the continuum subdomain) , 

+ [Q,A Pij ,A eij ] . 

5. Remove all particles from the continuum subdomain. 

6. Linearly interpolate the continuum state to the present time, 

(n ex - n P ) U H + U n H ext 



U 



H 



rip + 1 



1. Particle flux tracking 

As particles are advected during a particle time step, they may cross from a particle to a 
continuum cell and vice versa, and we need to keep track of the resulting fluxes. Note that a 
particle may cross up to d cell interfaces near corners, where d is the spatial dimension, and may 
even recross the same interface twice near hard-wall boundaries. Therefore, ray tracing is the most 
simple and reliable way to account for all particle fluxes correctly. For the majority of the particles 
in the interior, far from corners and hard walls, the usual quick DSMC update will, however, be 
sufficient. 

Algorithm 3: Move the fluid particle i from time t to time t + Atp and determine whether it 
crosses the coupling interface. We will assume that during a particle time step no particle can 
move more than one macro cell length along each dimension (in practice particles typically move 

only a fraction of a micro cell) . 

1. Store information on the macro cell c id to which the particle belongs, and then tentatively update 
the position of the particle <— r< + ViAtp, and find the tentative macro cell a and micro cell bi to 
which the particle moves, taking into account periodic boundary conditions. 

2. If c id and a are near a boundary, go to step [3] If c id = c%, or if c id and Cj are both continuum or 
are both particle cells and at least one of them is not at a corner, accept the new particle position 
and go to step [4] 

3. Undo the tentative particle update, <— r, — ViAtp, and then ray trace the path of the particle 
during this time step from one macro cell-cell interface I c to the next, accounting for boundary 
conditions (e.g., wrapping around periodic boundaries and colliding the particle with any hard walls 
it encounters |30j). Every time the particle crosses from a particle to a continuum cell, update the 
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particle flux at the cell interface I c , & P <— & p + [m, mvi, mvf /2] , similarly, if the particle crosses 
from a continuum to a particle cell, tfrp 1 ^ <— ffrp^ — \m, mvi, mvf /2l . 

4. If the new particle micro cell bi is neither in the particle subdomain nor in the reservoir region, remove 
the particle from the system. 



2. Inserting reservoir particles 

At every particle time step, reservoir particles need to be inserted into the reservoir region or 
unfilled continuum cells. These particles may later enter the particle subdomain or they may be 
discarded, while the trial particles in unfilled cells will be retained unless they leave the particle 
subdomain (see Algorithm [3]) . When inserting particles in an unfilled macro cell, it is important 
to maintain strict momentum and energy conservation by ensuring that the inserted particles have 
exactly the same total momentum and energy as the previous continuum values. This avoids global 
drifts of momentum and energy, which will be important in several of the examples we present in 



Section IV It is not possible to ensure strict mass conservation because of quantization effects, but 
by using smart (randomized) rounding one can avoid any spurious drifts in the mass. 

The velocity distribution for the trial particles should, at first approximation, be chosen from a 
Maxwell-Boltzmann distribution. However, it is well-known from kinetic theory that the presence 
of shear and temperature gradients skews the distribution, specifically, to first order in the gradients 
the Chapman- Enskog distribution is obtained |4"T] . It is important to note that only the kinetic 
contribution to the viscosity enters in the Chapman-Enskog distribution and not the full viscosity 
which also includes a collisional viscosity for dense gases jl6]. Previous work on deterministic 
DSMC hybrids has, as expected, found that using the Chapman-Enskog distribution improves 
the accuracy of the hybrids [20, 41 . However, in the presence of transient fluctuations and the 
associated transient gradients it becomes less clear what the appropriate distribution to use is, as 



we observe numerically in Section IV A 



Certainly at pure equilibrium we know that the Maxwell-Boltzmann distribution is correct de- 
spite the presence of fluctuating gradients. At the same time, we know that the CE distribution 
ought to be used in the presence of constant macroscopic gradients, as it is required in order to 
obtain the correct kinetic contribution to the viscous stress tensor [JH]. The inability to estimate 
time-dependent mean gradients from just a single fluctuating realization forces the use of instanta- 
neous gradients, which are unreliable due to the statistical uncertainty and can become unphysically 
large when Nq is small (say Nq < 50). In some cases the background macroscopic gradients may 



20 



be known a priori (for example, Couette or Fourier flows) and they can be used without trying to 
numerically estimate them, otherwise, they can be assumed to be zero or numerically estimated by 
performing some spatio-temporal averaging |21| . Note that here we assume that the particles are 
uniformly distributed (Poisson spatial distribution) as in an ideal gas, and we do not try to take 
the density gradient into account (but see the procedure described in the Appendix of Ref . [21] ) . 

Algorithm 4: Insert trial particles in the reservoir or unfilled macro cell c with centroid r c , 
taking into account the target density p c , velocity v c and temperature T c in cell c, calculated 
from the conserved cell state U c . If using the Chapman-Enskog distribution, also take into 
account the estimated local shear rate V c v and local temperature gradient V C T. 

1. Build a list C rc of the N rc micro cells contained in the macro cell c that need to be filled with particles. 
This typically excludes micro cells that are partially covered by impermeable beads or boundaries so 
as to avoid generating overlaps. 

2. Determine the total number N p of trial particles to insert into the reservoir portion of c by sampling 
from the binomial distribution 



where N p = \_p c V c /m\ is the total expected number of particles in c, V c is the volume of c, and 
p = N rc /N sc , where iV sc is the number of micro subcells per macro cell. For sufficiently large N p this 
can be well- approximated by a Gaussian distribution, which can be sampled faster. 

3. For each of the N p trial particles to be inserted, do: 

(a) Choose a micro cell b uniformly from the N rc micro cells in the list C rc . 

(b) Generate a random particle position r.; <— r c + r re i uniformly inside micro cell b. 

(c) Generate a random relative velocity for the particle v re i from the Maxwell-Boltzmann or 
Chapman-Enskog distribution [47] . and set the particle velocity by taking into account the de- 
sired continuum state in cell c and, if available, its estimated gradient, Vi <— v c + (V c u) r re i + 



(d) If the cell c is an unfilled macro cell, keep track of the total momentum P and energy E of the 
N p trial particles, to be adjusted in Step H] for conservation. 

4. If cell c was unfilled and N p > 1, then correct the particle velocities to match the desired total 
momentum P c — p c V c and energy E c = V c e c inside macro cell c, thus ensuring exact conservation: 
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(a) Calculate the scaling factor 

2 _ E c - P?/(2mN p ) 



a 



E - P 2 /(2mN p ) 
and velocity shift Av = (P — P c ) / (amN p ). 

(b) Scale and shift the velocity for every trial particle i, Vi <— a — Av). 



IV. RESULTS 

In this section we provide extensive tests of the hybrid scheme, in both equilibrium and non- 
equilibrium situations, and in both two and three dimensions. Our goal is to access how well the 
hybrid method can reproduce results obtained with a pure particle method, which we consider the 
gold standard. 

We have implemented our hybrid method in a code that can handle both two and three dimen- 
sional systems. Of course, one can study one- and two-dimensional flows with the three dimensional 
code by using periodic boundaries along the remaining dimensions. We refer to this as quasi one- 
or two-dimensional simulations. At the same time, both (I-)DSMC and the LLNS equations have 
a truly two-dimensional formulation, which we have also implemented for testing purposes. Trans- 
port coefficients for two-dimensional particle models formally diverge in the infinite-time limit (see 
the discussion in Ref. jl9]) and it is not obvious that the Navier-Stokes equations are a proper 
coarse-graining of the microscopic dynamics. However, this divergence is very slow (logarithmic) 
and it will be mollified (bounded) by finite system size, and we will therefore not need to con- 
cern ourselves with these issues. In the first three examples, we use the three-dimensional particle 
and continuum codes, and use the two-dimensional code only for the adiabatic piston example for 
computational reasons. 

We have also implemented continuum solvers for both the full non-linear and the linearized 
LLNS equations. As discussed in more detail in Ref. [9], the nonlinear LLNS equations are 
mathematically ill-defined and this can lead to breakdown in the numerical solution such as negative 
densities or temperatures. At the same time, the linearized equations are not able to describe a 
wide range of physical phenomena such as the effect of fluctuations on the mean flow; they also omit 
a number of terms of order Nq 1 , where Nq is the average number of particles per continuum cell 
(e.g., the center-of-mass kinetic energy). If the number of particles per continuum cell is sufficiently 
large (in our experience, Nq > 75) the fluctuations are small and the difference between the linear 
and nonlinear hydrodynamic solvers is very small, and we prefer to use the nonlinear solver. We 
will use Nq ~ 100 in our hybrid simulations. 
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In our implementation, the continuum solver can either be the more accurate third-order Runge- 
Kutta (RK3) temporal integrator developed in Ref. [9] or a more classical stochastic MacCormack 
integrator [BJ. The analysis in Ref. [§] shows that obtaining reasonably-accurate equilibrium 
fluctuations with predictor-corrector methods for the diffusive fluxes, as used in the MacCormack 
scheme, requires using a continuum time step Ate that is a fraction of the CFL stability limit 
AtcFL- In the simulations we present here we have typically used a time step At ~ 0.2AtcFL, 
which is typically still about 5 times larger than the particle time step Atp, and we have found 
little impact of the exact value of Ate- 

The hybrid method requires estimates of the transport coefficients of the particle fluid, notably 
the viscosity and the thermal conductivity. For traditional DSMC at low densities there are rather 
accurate theoretical estimates of the viscosity and thermal conductivity [50} [51], however, it is 
nontrivial to obtain reasonably accurate theoretical values at the higher densities we use in I- 
DSMC because of the importance of multi-particle correlations. We therefore estimate the transport 
coefficients numerically. For this purpose we simulate a system with periodic boundaries along two 
of the directions and isothermal stick wall boundaries along the other direction. To estimate the 
viscosity we apply a shear flow by moving one of the wall boundaries at constant speed inducing 
a Couette flow with an approximately constant shear gradient Vi; = ^(Vi> + Vi> r ). We then 
calculate the steady-state stress tensor cr by averaging over all particles i and colliding pairs of 
particles ij over a long time interval At, 



cr = cr k + <r c = m (vi (g> Vi) + 



At 



where <Xfc = PI + 2r}k^v is the kinetic contribution giving the kinetic viscosity and cr c = 2r/ c Vi> 
is the collisional contribution giving the collisional viscosity rj c , r\ = % + rj c . We exclude particles 
that are close to the wall boundaries when calculating these averages to minimize finite-size effects. 
Similarly, for the thermal conductivity we apply a small constant temperature gradient VT by 
setting the two walls at different temperatures, and we also impose the required density gradient 
to maintain mechanical equilibrium (constant pressure). We then calculate the steady-state heat 
flux vector £ = //VT, 



.27' At 

from which we obtain the kinetic and collisional contributions to the thermal conductivity. There 
are alternative methods that one can use to calculate the transport coefficients, using both equilib- 
rium and non-equilibrium settings, however, we have found the above method to be most accurate 
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for a given computational effort if only moderate accuracy is desired. Results from different non- 
equilibrium methods are found to be within 5 — 10% of each other. 

A. Mismatch at the interface 

Previous work has studied a hybrid scheme very similar to the one described here for several 
quasi one-dimensional situations [21_. Reference [21] first studied a pure equilibrium situation in 
which one part of a periodic domain was covered by a particle subdomain, and found that the 
stochastic hybrid scheme was able to reproduce the spatio-temporal correlations in equilibrium 
fluctuations very well, although some mismatch at the particle/continuum interface was found. 
Here we explore this mismatch more carefully by studying a quasi one-dimensional periodic system 
where the middle portion of the domain is filled with particles and the rest is continuum. In each 
macro cell, we compute the average density, temperature and velocity and their variances with 
high accuracy. 

We first calculate the mean conserved quantities (mean density, momentum density, and energy 
density) in each macro cell, and then calculate the mean velocity and temperature from those, for 
example, (v) = (j) / (p) instead of averaging the instantaneous velocities, (v) = (j/p). As shown 
in Ref. [52] , the approach we use leads to an unbiased estimate of the mean, while the latter has a 
bias when there are correlations between the fluctuations of the different hydrodynamic variables. 
The variances are estimated from the instantaneous values, e.g., (5v 2 ^ = (^(j / p) 2 ^ — (j/p) 2 - It is 
possible to construct unbiased estimates for the variances as well [52], however, this is somewhat 
more involved and the bias is rather small compared to the artifacts we are focusing on. 

In Fig. [3] we show the means and variances along the length of the system, normalized by 
the expected values [6]. For velocity, the mean is zero to within statistical uncertainty in both 
the particle and continuum subdomains, and we do not show it in the figure. However, a small 



mismatch is clearly seen in Fig. 3(a) between the density and temperature in the particle and 
continuum subdomains. The mismatch is such that the pressure (p) = {p) R (T) is constant across 
the interface, that is, the particle and continuum subdomains are in mechanical equilibrium but not 
in true thermodynamic equilibrium. In Appendix [A] we show that this kind of mismatch is expected 
because the average particle fluxes coming from the reservoir particles inserted in the continuum 
subdomains have a bias of order Nq 1 , where Nq is the average number of particles per macroscopic 
cell. Because our coupling matches both the state and the fluxes across the interface this bias 
makes it impossible for the particle and continuum to reach true thermodynamic equilibrium. The 
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theory in Appendix |a| suggests that the size of the mismatch is of order Nq 1 , consistent with 
the results in Fig. |3(a)| however, the crude theoretical estimates do not actually give the steady 
state since they assume equilibrium to begin with. That the cells near the interface are not in 
equilibrium is reflected in the variances of the hydrodynamic variables, which have a spike near the 
particle-continuum interface, as shown in Fig. |3(b)| We have observed that the relative magnitude 
of this spike does not depend on iVo- 

The cause of the mismatch is the fact that we use the instantaneous values of the local density, 
velocity and temperature when generating the velocities for the reservoir particles. This is necessary 
because we cannot in general obtain estimates of the time-dependent mean values from running a 
single realization, and are forced to use the instantaneous values. One can use some sort of spatio- 
temporal averaging to obtain estimates of the local means, however, this introduces additional 
time and length scales into the algorithm that do not have an obvious physical interpretation. For 
steady-state problems it may be possible to use running means to avoid the mismatch we observe, 
however, this is not possible in a general dynamic context. Deeper theoretical understanding of 
the connection between the microscopic dynamics and the LLNS equations is necessary to design 
a more consistent approach. 

Another complex issue that arises in the fluctuating hybrid is whether to use the Maxwell- 
Boltzmann (MB) or the Chapman-Enskog (CE) distributions when generating the velocities of 



the reservoir particles (see discussion in Section HID 2). In Fig. 3(a) we compare the size of the 
mismatch at the interface when the MB and CE distributions are used. For the CE distribution 
we obtained a local estimate of the gradient using simple centered differences of the instantaneous 
hydrodynamic variables. As seen in the figure, there is a greater discrepancy when the CE dis- 
tribution is used, especially for the variances. We will therefore adopt a compromise in which we 
use the MB distribution for all calculations reported here. In cases when there is a macroscopic 
background gradient that is specified a priori (e.g., shear flows) we can use that gradient in the 
CE distribution. 

Note that we have performed numerous additional quasi-one dimensional tests of the coupling 
that we do not describe here for brevity. For example, we have tested the matching of the shear 
stress tensor in a shear flow parallel to the particle-continuum interface by verifying that a linear 
velocity profile is obtained without any slope discontinuity at the interface (see Appendix D.2 in 
Ref. [S]). We have also reproduced the results reported in Ref. [21J, such as the presence of 
unphysical long-range correlations in the fluctuations in the particle region when the deterministic 
instead of the stochastic hybrid is used. 
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Figure 3: Normalized means and variances of the hydrodynamic variables in each of the 46 macroscopic 
cells of a quasi one-dimensional periodic system where the middle portion (20 continuum cells, in-between 
dashed vertical lines) is the particle subdomain. The velocities of the reservoir particles are either samples 
from a Maxwell-Boltzmann (MB, lines only) or a Chapman-Enskog distribution (CE, lines and symbols), 
(a) Unbiased [55] estimates of the mean density and temperature (we averaged about 2.5 • 10 6 samples taken 
every macro step) for two different sizes of the continuum cells, ones containing A^o = 120 particles on 
average, and ones containing No — 480 particles. Only the MB results are shown for the larger cells for 
clarity, with similar results observed for CE. |(b)| Estimates of the variances of the hydrodynamic variables 
for No — 120 particles per continuum cell (also an average over 2.5 • 10 6 samples). 
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B. Dynamic Structure Factors 

The hydrodynamics of the spontaneous thermal fluctuations in the I-DSMC fluid is expected 
to be described by the Landau-Lifshitz Navier-Stokes (LLNS) equations for the fluctuating field 
U = (po + Sp,6v,To + 5T) linearized around a reference equilibrium state Uq = (pq,vq = 0,Tq) 
[35j E3] • By solving these equations in the Fourier wavevector-frequency (k,ui) domain for 
U(k,uj) = (5p,5v,5T) and performing an ensemble average over the fluctuating stresses we can 
obtain the equilibrium (stationary) spatio-temporal correlations (covariance) of the fluctuating 
fields. We express these correlations in terms of the 3x3 symmetric positive-definite hydrodynamic 
structure factor matrix Sn(k,uj) = (uU \ [9 , which is essentially the spatio-temporal spectrum 
of the fluctuating fields. Integrating Sn(k, uj) over frequencies gives the hydrostatic structure factor 
matrix Sn(k), which turns out to be diagonal since in any given snapshot the hydrodynamic 
variables are uncorrelated at equilibrium. 

We non-dimensionalize Sn(k, uj) so that Sn{k) is the identity matrix. For example, the density- 
density correlations are given by the dimensionless structure factor 

S p (k,uj) = (poCQ 2 k B T ) 1 (p(k,uj)p*(k,uj)) , 

and we express the spatio-temporal cross-correlation between density and velocity through the 
dimensionless structure factor 

S PtV (k,<J) = (p Co 2 ^7bp (po'^Top {p(k,u)v*{k,u)} , 

where Cg = ksTQ/m is the isothermal speed of sound. Reference [21] demonstrated that a hybrid 
method very similar to the one we described here correctly reproduces the density-density time- 
correlation function S p (k, t) for large wavelengths in a quasi one-dimensional periodic system. The 
density-density dynamic structure factor S p (k, uj) is often the only one considered because it is 
accessible experimentally via light scattering measurements and thus most familiar [53]. The full 
dynamic structure factor matrix Sn(k, uj) is a more complete measure of the spatio-temporal evo- 
lution of the thermal fluctuations and includes both sound (hyperbolic) and dissipative (diffusive) 
effects. It is therefore important to show that the hybrid scheme correctly reproduces Sji(k,uj) as 
compared to a purely particle simulation and demonstrate that the hybrid is capable of capturing 
the propagation of spontaneous thermal fluctuations across the particle-continuum interface. 
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1. Bulk Dynamic Structure Factor 

For a bulk fluid, i.e., for periodic boundary conditions, it is well-known |53 j that the density- 
density component S p (k,u) and the temperature-temperature component Sr(fe>w) of Sn(k,uj) 
exhibit three peaks for a given wavevector k. There is one central Rayleigh peak at u = which 
comes from entropic fluctuations. There is also two symmetric Brillouin peaks at uo ~ c s k, where 
c s is the adiabatic speed of sound, which come from the isoentropic propagation of sound waves 
induced by the fluctuations. For the components of the velocity parallel to the wavevector the 
dynamic structure factors S v ,, (k, to) exhibit all three peaks, while for the component perpendicular 
to the wavevector S v± (k,uj) lacks the central peak. 

Figure [4] shows selected dynamic structure factors for a quasi two-dimensional system with 
periodic boundary conditions. The simulation box is composed of 10x10x1 macro cells, each cell 
containing about 120 particles. Finite-volume averages of the hydrodynamic conserved variables 
and the corresponding spatial Discrete Fourier Transforms (DFTs) were then calculated for each 
cell every 10 macro time steps and a temporal DFT was used to obtain discrete dynamic structure 
factors [9] for several wavenumbers. In the first set of particle runs, the whole domain was filled 
with particles and the macro cells were only used to sample hydrodynamic fields. In the second set 
of hybrid runs, the central portion of the simulation box was split along the x axes and designated 
as a particle subdomain, and the remaining two thirds of the domain were continuum. 

In Fig. [4] we show the results for a wavevector k that is neither parallel nor perpendicular to 
the particle-continuum interface so as to test the propagation of both perpendicular and tangen- 
tial fluctuations across the interface. The results show very little discrepancy between the pure 
particle and the hybrid runs, and they also conform to the theoretical predictions based on the 
LLNS equations. Perfect agreement is not expected because the theory is for the spectrum of the 
continuum field while the numerical results are discrete spectra of cell averages of the field, a dis- 
tinction that becomes important when the wavelength is comparable to the cell size. Additionally, 
even purely continuum calculations do not reproduce the theory exactly because of spatio-temporal 
discretization artifacts. 

2. Dynamic Structure Factors for Finite Systems 

The previous section discussed the "bulk" dynamic structure factors, as obtained by using peri- 
odic boundary conditions. For non-periodic (i.e., finite) systems equilibrium statistical mechanics 
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Figure 4: Discrete dynamic structure factors for waveindices q x — k x L x / (2-k) = 1, q y = 2 and q z — for a 
quasi two-dimensional system with lengths L x = L y = 2, L z = 0.2, split into a grid of 10 x 10 x 1 macro cells. 
Each macro cell contains about 120 I-DSMC particles of diameter D — 0.04 (density <f> = 0.5, and collision 
frequency \ — 0.62). An average over 5 runs each containing 10 4 temporal snapshots is employed. We 
perform pure particle runs and also hybrid runs in which only a strip 4 macro cells along the x axes was filled 
with particles and the rest handled with a continuum solver. The different hydrodynamic pairs of variables 
are shown with different colors, using a solid line for the result from the hybrid runs, symbols for the results 
of the pure particle runs, and a dashed line for the theoretical predictions based on the linearized LLNS 
equations (solved using the computer algebra system Maple). (Top) The diagonal components S p (k,ui), 
S Vx (k,uj), S Vy (k,uj) and Sr(fe,w). (Bottom) The off-diagonal components (cross-correlations) S PiV9 (k,u>), 
S VxtVy (k,oj), S Pt T(k,u) and S Vx>T (k,w). 
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requires that the static structure factor be oblivious to the presence of walls. However, the dynamic 
structure factors exhibit additional peaks due to the reflections of sound waves from the bound- 
aries. At a hard- wall boundary surface d£l with normal vector n either Dirichlet or von Neumann 
boundary conditions need to be imposed on the components of the velocity and the temperature 
(the boundary condition for density follows from these two). Two particularly common types of 
boundaries are: 

Thermal walls for which a stick condition is imposed on the velocity, vqq = 0, and the temper- 
ature is fixed, Tqq = Tq. 

Adiabatic walls for which a slip condition is imposed on the velocity, n ■ Vvn = and v± = 0, 
and there is no heat conduction through the wall, n ■ VT = 0. 

In particle simulations, these boundary conditions are imposed by employing standard rules for 
particle reflection at the boundaries [30 . We describe the corresponding handling in continuum 
simulations in Appendix [Cj In Appendix [B] we derive the form of the additional peaks in the 
dynamic structure factor for adiabatic walls by solving the linearized LLNS equations with the 
appropriate conditions. 

In Figure [5] we show dynamic structure factors for a quasi one-dimensional system bounded by 
adiabatic walls. As for the bulk (periodic) case in the previous section, we perform both purely 
particle runs and also hybrid runs in which the middle third of the domain is designated as a 
particle subdomain. Additional peaks due to the reflections of sound waves from the boundaries 
are clearly visible and correctly predicted by the LLNS equations and also accurately reproduced 
by both the purely continuum solution (not shown) and the hybrid. Similar agreement (not shown) 
is obtained between the particle, continuum and the hybrid runs for thermal walls. These results 
show that the hybrid is capable of capturing the dynamics of the fluctuations even in the presence 
of boundaries. Note that when the deterministic hybrid scheme is used one obtains essentially the 
correct shape of the peaks in the structure factor (not shown), however, the magnitude is smaller 
(by a factor of 2.5 for the example in Fig. [5| than the correct value due to the reduced level of 
fluctuations. 

C. Bead VACF 

As an illustration of the correct hydrodynamic behavior of the hybrid algorithm, we study 
the velocity autocorrelation function (VACF) C(t) = (v x (0)v x (t)) for a large neutrally-buoyant 
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Figure 5: Discrete dynamic structure factors for a quasi one-dimensional system with length L = 7.2 
(corresponding to 36 continuum cells) bounded by adiabatic walls, for waveindex q = 2 and wavevector 
k = 2irq/L » 1.75. The I-DSMC fluid parameters are as in Fig. [4] We perform pure particle runs and 
also hybrid runs in which the middle third of the domain is filled with particles and the rest handled with 
a continuum solver. The results from purely particle runs are shown with symbols, while the results from 
the hybrid are shown with a solid line. The theoretical predictions based on the linearized LLNS equations 
and the equations in Appendix [B] (solved using the computer algebra system Maple) are shown with a 
dashed line for adiabatic boundaries and dotted line for periodic boundaries. Since the magnitude of the 
Brillouin peaks shrinks to one half the bulk value in the presence of adiabatic walls, we also show the result 
with periodic boundaries scaled by 1/2 (dashed-dotted line). (Top) Dynamic structure factor for density, 
Sp(k,uj), showing the Raylcigh peak and the multiple Brillouin peaks. (Bottom) Dynamic structure factor 
for the component of velocity perpendicular to the wall, S v± (k,u)), which lacks the Rayleigh peak. The 
corresponding correlations for either of the parallel velocity components, S v ,.(k,u>), have only a Rayleigh 
peak, shown in the inset on a log-log scale. 
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impermeable bead of mass M and radius R diffusing through a dense Maxwell I-DSMC stochastic 
fluid [H] of particles with mass m <C M and collision diameter D <C R and density (volume 
fraction) (j) [mass density p = 6m(j>/(nD 3 )]. The VACF probl em is relevant to the modeling 
of polymer chains or (nano) colloids in solution (i.e., complex fluids), in particular, the integral 
of the VACF determines the diffusion coefficient which is an important macroscopic quantity. 
Furthermore, the very first MD studies of the VACF for fluid molecules led to the discovery of a 
long power-law tail in C{t) [53] which has since become a standard test for hydrodynamic behavior 
of methods for complex fluids [55j ESI ES EE1 E3 [60j EE] • 

The fluctuation-dissipation principle [62 points out that C(t) is exactly the decaying speed of a 
bead that initially has a unit speed, if only viscous dissipation was present without fluctuations, and 
the equipartition principle tells us that C(0) = {v%) = kT/2M. Using these two observations and 
assuming that the dissipation is well-described by a continuum approximation with stick boundary 
conditions on a sphere of radius Rh, C{t) has been calculated from the linearized (compressible) 
Navier-Stokes (NS) equations |63| 151] . The results are analytically complex even in the Laplace 
domain, however, at short times an inviscid compressible approximation applies. At large times 
the compressibility does not play a role and the incompressible NS equations can be used to predict 
the long-time tail |64[ [65] . At short times, t < t c = 2Rjj/c s , the major effect of compressibility 
is that sound waves generated by the motion of the suspended particle carry away a fraction of 
the momentum with the sound speed c s , so that the VACF quickly decays from its initial value 
C(0) = k B T/M to C(t c ) « k B T/M e ff, where M eff = M + 2irR 3 p/3 [63]. At long times, t > t visc = 
ApR 2 H /?>ri, the VACF decays as with an asymptotic power-law tail (fcgT '/M)(8\ /r 3Tr)~ 1 (t/t v i sc )^ 3 ^ 2 , 
in disagreement with predictions based on the Langevin equation (Brownian dynamics), C{t) = 
(k B T/M) exp (-6TrR H r]t/M). 

We performed purely particle simulations of a diffusing bead in various I-DSMC fluids in Refs. 
|18[ [6"T] . In purely particle methods the length of the runs necessary to achieve sufficient accuracy 
in the region of the hydrodynamic tail is often prohibitively large for beads much larger than the 
fluid particles themselves. It is necessary to use hybrid methods and limit the particle region to the 
vicinity of the bead in order to achieve a sufficient separation of the molecular, sonic, viscous, and 
diffusive time scales and study sufficiently large box sizes over sufficiently long times. In the results 
we report here we have used an impermeable hard bead for easier comparison with existing theory; 
similar results are obtained using permeable beads. The interaction between the I-DSMC fluid 
particles and the bead is treated as if the fluid particles are hard spheres of diameter D s , chosen 
to be somewhat smaller than their interaction diameter with other fluid particles (specifically, 
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we use D s = D/4) for computational efficiency reasons, using an event-driven algorithm [30 . 
Upon collision with the bead the relative velocity of the fluid particle is reversed in order to 
provide a no-slip condition at the surface of the suspended sphere [6U\ 158] (slip boundaries give 
qualitatively identical results). We have estimated the effective (hydrodynamic) colloid radius 
Rh from numerical measurements of the Stokes friction force F = —Qt^Rh^v and found it to be 
somewhat larger than the hard-core collision radius R + D s /2, but for the calculations below we 
use Rh = R + D s /2. Since we used periodic boundary conditions with a box of length L, there are 
artifacts in C{t) after about the time at which sound waves generated by its periodic images reach 
the particle, = L/c s . The averaging procedure that we used in order to eliminate some of the 
noise from the tail does not properly resolve these sound effects 2 , even though they are visible in 
the results shown below. Furthermore, there are strong finite-size effects that manifest themselves 
as a rapid decay of the VACF for times longer than the viscous time-scale pL 2 /rj |59|. 

For the hybrid calculations, we localize the particle subdomain to the continuum cells that 
overlap or are close to the moving bead. The location of the particle subdomain is updated 
periodically as the bead moves. The algorithm that we use tries to fit the particle subdomain as 
closely around the bead but ensuring that there are a certain number of micro cells in-between the 
surface of the bead and the particle-continuum interface. The exact shape of the particle subdomain 
thus changes as the bead moves and the number of particles employed by the hybrid fluctuates, 
especially when the bead is small compared to the continuum cells. Obtaining reasonably-accurate 
results for the VACF at long times requires very long runs. We found that it is crucial to strictly 
conserve momentum in the hybrid when unfilled continuum cells are transformed into particle 
cells. Otherwise very slow drifts in the momentum of the system appear due to the use of periodic 
boundary conditions, and this drift changes the tail of C(t), especially for massive beads where 
the typical bead speed is already small compared to the typical fluid particle speed. We used the 
MacCormack solver and the linearized formulation of the LLNS equations for these simulations, 
however, similar results are obtained with the non-linear Runge-Kutta solver as long as the macro 
cells are sufficiently large. As discussed earlier, the use of the linearized formulation makes it 
possible to reduce the size of the continuum cells without introducing numerical problems due to 

2 To calculate the VACF, we first calculated the average mean-square displacement Ar(t) by averaging over a long 
run and numerically differentiated this to obtain a time-dependent diffusion coefficient D(t). We then smoothed 
D(t) by fitting a quadratic polynomial over short-time intervals spaced on a logarithmic scale (so that more points 
were averaged over in the tail), and obtained the VACF by differentiating the smoothed D(t). This procedure 
produces well-resolved tails, however, it obscures features over short time scales compared to the smoothing interval. 
A more direct velocity autocorrelation calculation was used at very short times. 
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the non-linearity of the equations. 




(a) 



(b) 



Figure 6: Normalized velocity autocorrelation function (VACF) C(t) / (ksT/M) for a neutrally-buoyant hard 
spherical bead of mass M suspended in a fluid of I-DSMC particles of diameter D, for two different bead 
sizes, a small bead of radius R = 1.2 5D [(a)] and a large bead of radius R = 6.25£) [(b)] A log-log scale is used 
to emphasize the long-time power law tail and the time is normalized by the viscous time t V i SC , so that the 
results should be approximately independent of the actual bead radius. The inset shows the initial decay of 
the VACF on a semi-log scale, where the time is normalized by the sonic time scale t c . Periodic boundary 
conditions with a cubic cell of length L are employed, and the sound crossing time ti is indicated. Results 
from purely particle simulations are shown with a dashed-dotted line, and the incompressible hydrodynamic 
theory is shown in with a dotted line. Results from hybrid runs are also shown with a solid line for the 
stochastic hybrid and dashed line for the deterministic hybrid, for both the same box size as the particle 
runs (red) and a larger simulation box (green). 



In the calculations reported in Fig. [6j the I-DSMC fluid has a density (volume fraction) <fi = 0.5 
and collision frequency prefactor x = 0.62. The adiabatic sound speed is c s = \Jhk^T j (3m) and 
viscosity is t] = r)D~ 2 \J raksT , where we measured r\ ~ 0.75. Note that in atomistic time units 
to = D^/m/ksT the viscous time scale is t v i sc /to 6(()(Rh / D) 2 / (3irfj) a* 0A(Rh /D) 2 . 

As a first test case, in Fig. |6(a)| we compare against the particle data from Ref. |61j . for which 
the size of the bead is R = 1.25D, M = 7.81m, and the simulation box is L = 1 = 25D, which 
corresponds to 24 3 micro cells and about N ~ 1.5 • 10 4 particles. The hybrid runs used macro 
cells each composed of 4 3 micro cells, which corresponds to about A^o = 80 particles per cell, and 
the size of the particle subdomain fluctuated between about 3 • 10 3 and 6 • 10 3 particles due to the 



change of the location of the bead relative to the continuum grid. The particle result in Fig. 6(a) 
is the average over 10 runs, each of length T/t v i sc ~ 2 • 10 5 , while the hybrid results are from a 
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single run of length T/t v i sc ~ 7.5 • 10 5 . It is seen in the figure that both the deterministic and the 
fluctuating hybrid reproduce the particle results closely, with a small but visible difference at long 
times where the deterministic hybrid under-predicts the magnitude of the tail in the VACF. We 
also show results from a hybrid run with a twice larger simulation box, L = 2, which marginally 
increases the computational effort in the hybrid runs, but would increase the length of purely 
particle runs by an order of magnitude. The hydrodynamic tail becomes pronounced and closer to 
the theoretical prediction for an infinite system, as expected. We have observed (not shown) that 
using small continuum cells composed of only 3 3 micro cells, which corresponds to about Nq = 35 
particles per cell, leads to an over-prediction of the magnitude of the VACF at short times, that 
is, to an excess kinetic energy for the bead by as much as 20%, depending on the exact parameters 
used. 

As a second, more difficult test, in Fig. |6(b)| we report results from particle simulations for 
a much larger bead, R = 6.25D, M = 976m, and the simulation box is L = 2 = 50D, which 
corresponds to N ~ 1.2 ■ 10 5 particles. We have performed a variety of hybrid runs and in Fig. 



6(b) we report results from runs with macro cells each 3 micro cells, as well as results for a larger 



simulation box, L = 3, and macro cells each composed of 4 3 micro cells. The particle results are 
the average over 5 runs each of length T/t v i sc ~ 2.5 • 10 3 , while the hybrid results are from a single 
run of length T/t v i sc ~ 7.5 • 10 3 . The hybrid runs had a particle subdomain containing about 2 • 10 
particles. We observed little impact of the size of the continuum cell size or the size of the particle 
subdomain. The results in Fig. |6(b)| show that the stochastic hybrid correctly reproduces the tail 
in the VACF, while it slightly over-estimates the VACF at short times. The deterministic hybrid, 
on the other hand, strongly under-estimates the magnitude of C(t) at both short and long times. 
It is particularly striking that the deterministic hybrid fails to reproduce the magnitude of the long 
time tail (and thus the diffusion coefficient), vividly demonstrating the importance of including 
fluctuations in the continuum domain. 

Computing the VACF for a diffusing bead has become a standard test for micro- and nano-scale 
fluid-structure coupling methods and has been performed for a suspended bead in a wide range of 
particle and (semi-)continuum compressible and incompressible fluids [SH |S3 [Ml E3 EH1 EH1 USB EI] • 
However, these tests often do not correctly test for all of the components necessary to match the 
VACF over all relevant timescales: equipartition, a fluctuation-dissipation relation, and hydrody- 
namics. Purely continuum fluid methods allow for using a much larger time step than particle 
(and thus hybrid) methods, especially if an incompressible formulation is directly coupled to the 
equations of motion of the suspended bead [53 |59l 60 J. When fluctuations are not included in 
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continuum methods, the VACF is often obtained by considering the deterministic decay of the ve- 
locity of a bead. This however assumes a priori that proper thermodynamic equilibrium exists with 
the correct fluctuation-dissipation relation, without actually testing this explicitly. Alternatively, 
coupling methods between a continuum fluid and a suspended particle often have some arbitrary 
coupling parameters that are tuned to reproduce the desired diffusion coefficient without producing 
a physically-consistent VACF, especially at short times In particular, incompressible formu- 
lations cannot reproduce the initial value or the decay of the VACF and should instead aim to 
produce an average kinetic energy of the bead of ksT rather than the statistical mechanics result 
of 2>UbT/2 [65]. It is therefore important to more carefully access the ability of other methods in 
the literature to correctly reproduce the full VACF for a truly equilibrium simulation with a large 
bead. 

D. Adiabatic Piston 

The problem of how thermodynamic equilibrium is established has a long history in statistical 
mechanics [HZ] . The adiabatic piston problem is one of the examples used to study the fluctuation- 
driven relaxation toward equilibrium |68[ 169) \7U\ ITT] that is simple enough to be amenable to 
theoretical analysis but also sufficiently realistic to be relevant to important problems in nano- 
science such as Brownian motors |24l 172] . We study the following formulation of the adiabatic 
piston problem. A long quasi one-dimensional box with adiabatic walls is divided in two halves with 
a thermally-insulating piston of mass M 3> m that can move along the length of the box without 
friction. Each of the two halves of the box is filled with a fluid that is, initially, at a different 
temperature T and density p, here assumed to follow the ideal equation of state P = pk^T jm. 
If the macroscopic pressures in the two halves are different, plTl PrTr, then the pressure 
difference will push the piston to perform macroscopic oscillations with a period that can be 
estimated by assuming that each half undergoes an adiabatic transformation {PV" 1 = const.). 
These oscillations are damped by viscous friction and lead to the piston essentially coming to rest in 
a state of mechanical equilibrium, PlTl = PrTr. This stage of the relaxation from non-equilibrium 
to mechanical equilibrium has been shown to be well-described by deterministic hydrodynamics 

The state of mechanical equilibrium is however not a state of true thermodynamic equilibrium, 
which also requires equality of the temperatures on the two sides of the piston. Reaching full 
equilibrium requires heat transfer through the piston, but the piston is adiabatic and does not 
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conduct heat. In classical deterministic hydrodynamics the piston would just stand still and never 
reach full equilibrium. It has been realized long ago that heat is slowly transferred through the 
mechanical asymmetric fluctuations of the piston due to its thermal motion, until the temperatures 
on both sides of the piston equilibrate and the fluctuations become symmetric. This equilibration 
takes place through a single degree of freedom (the piston position) coupling the two large reservoirs, 
and it would be astronomically slow in a laboratory setting. While various Langevin or kinetic 
theories have been developed for the effective heat conduction of the adiabatic piston (see Refs. 
|68|. 169] ITOT ITT] and references therein), there is no complete theoretical understanding of the 
effective heat conductivity, especially in dense fluids. Molecular dynamic simulations have been 
performed in the past |68| 169] [70] using hard-disk fluids, but the very long runs required to reach 
thermodynamic equilibrium for massive pistons have limited the size of the systems that could 
be studied. Furthermore, there have been no studies applying fluctuating hydrodynamics to this 
problem. 

Here we apply our hybrid method to the adiabatic piston problem in two dimensions, using a 
non-linear two-dimensional implementation of the Runge-Kutta integrator described in Ref. [9] as 
the continuum solver. The choice of two dimensions is for purely computational reasons. Firstly, 
the number of particles required to fill a box of sufficient size is much smaller thus allowing for long 
particle simulations. Secondly, in order to implement the piston in our particle scheme we reused 
the same mixed event-driven/time-driven handling [30] as we used for the VACF computations in 
Section IV C| Namely, we made a piston out of Nb small impermeable beads, connected together 
to form a barrier between the two box halves, as illustrated in Fig. [7] In two dimensions, by 
ensuring that two piston beads never separate by more than a given distance we can ensure that 
two I-DSMC particles on opposite sides of the piston cannot possibly collide and thus the piston 
will be insulating. We have studied two different types of pistons, a flexible piston where the 
the beads are tethered together to form a chain [30J that is stretched but where each individual 
bead can still move independently of the others, and a rigid piston that is obtained with a slight 
modification of the event loop to move all of the piston beads in unison. While at the macroscopic 
level the exact shape of the piston should not make a big difference, we have found that increasing 
the number of degrees of freedom of the piston from one to iVj, makes a significant difference 
in the thermal conductivity of the piston, and therefore, we will focus here on rigid pistons as 
in the traditional formulation. We use specular collisions of the fluid particles with the piston 
beads, although qualitatively identical (but not quantitatively identical) results are obtained using 
bounce-back collisions as well. 
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Figure 7: An illustration of the computational setup used for the adiabatic piston computations. Only the 
central portion of the box of aspect ratio 6 x 1 is shown. Left of the piston the gas is cold and dense; to the 
right it is hot and dilute. The piston beads (red disks) separate the box into two halves, and are surrounded 
on each side by a fluid of I-DSMC particles (smaller green disks), which is twice denser but also twice cooler 
in the left half than in the right half. The microscopic grid is shown with thinner light blue lines and the 
hydrodynamic grid is shown with thicker dark blue lines. The interface between the particle and continuum 
regions is highlighted with a thick red line. A snapshot of the values of the hydrodynamic variables in each 
continuum cell is shown using a large purple disk whose size is proportional to the density and its opaqueness 
is proportional to the temperature, and an arrow for the fluctuating velocities. 

The hybrid method setup for the adiabatic piston is illustrated in Fig. [7j We use a two- 
dimensional Maxwell I-DSMC particle fluid (0 = 1, x = 1) with collision diameter D = 0.1 (hard- 
sphere diameter D s = D/2) and a piston composed of N\, = 40 beads of diameter Df, = 0.0955. The 
particle subdomain is limited to a few continuum cells around the piston, which we keep at about 
two or more continuum cells on each side of the piston, so that the unreasonable hydrodynamic 
values in the cells that overlap the piston do not affect the continuum solver appreciably. Periodic 
boundary conditions are applied along the y dimension (parallel to the piston) with the width of 
the domain L y = 4 being 40 microscopic cells, while adiabatic walls were placed at the ends of 
the box whose total length L x = 24 was 240 microscopic cells. We have studied various sizes for 
the macroscopic cells, and report results for a quasi one-dimensional continuum grid in which each 
macro cell contains 4 x 40 micro cells, corresponding to about 200 particles per continuum cell. We 
also present results for a two-dimensional continuum grid where each macro cell contains 8 x 10 
micro cells. 

We have performed hybrid runs with both the deterministic and stochastic hybrids. In Fig. |8]we 
show the position of a massive piston of mass M = 4000m that started at a position x = 8 that is 
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Figure 8: Relaxation of a massive rigid piston (M/m = 4000) from a position x = 8 that is not in mechanical 
equilibrium. Through rapid damped oscillations mechanical equilibrium is established at position x ~ 7, 
after which a slow relaxation to true equilibrium is seen. The stochastic hybrid is able to match the particle 
data very well, to within the expected statistical difference from realization to realization. The deterministic 
hybrid, on the other hand, clearly under-predicts the rate of relaxation. The inset highlights the initial 
oscillations and shows that in the regime where fluctuations do not matter the deterministic and stochastic 
hybrids do not differ appreciably. 

not in mechanical equilibrium and thus performs rapid damped oscillations until it reaches a state of 
mechanical equilibrium at x ~ 7, from which it slowly relaxes toward true equilibrium. The results 
in the figure show that the stochastic hybrid reproduces the correct relaxation toward equilibrium 
while the deterministic hybrid severely under-predicted the rate of equilibration (effective heat 
conductivity), even though the initial mechanical stage of the relaxation is correctly captured by 
both hybrids, as expected. We have observed that the deterministic hybrid fails to give the correct 
answer whenever a rigid massive piston is used, M > 250m. For flexible pistons, we find that even 
for a large bead mass M;,(overall piston mass M = N^M^) both the deterministic and stochastic 
hybrids reproduce the purely particle results for the slow relaxation toward equilibrium. 

A more detailed comparison of the particle and hybrid results for a piston of mass M = 1000m 
that is initially in mechanical equilibrium at position x = 8 = L x /3 is shown in Fig. [9] The initial 
conditions were ksTh = 2/3, pl = 2/3 and /cbTr = 4/3, pl = 1/3, so that there is an equal mass 
on each side of the piston. At the true equilibrium state the piston remains close to the middle 
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Figure 9: Relaxation of a rigid piston of mass M/m = 1000 from an initial state of mechanical equilibrium 
[x = 8) to a state of thermodynamic equilibrium (x — 12). The inset emphasizes the initial exponential 
decay on a semi- log scale. The hybrid runs used a particle subdomain of width wp = 2 on each side of 
the piston and continuum cells that were composed of either 4 x 40 or 8 x 10 microscopic cells. For the 
deterministic hybrid the macro cell size makes little difference so we only show the 4 x 40 case. 



of the box, x eq = L/2 = 12, with equal density on each side. The results shown are averages over 
10 samples, but it should be emphasized that each run exhibits thermal oscillations of the piston 3 
position that are diminished by direct averaging because they have different (random) phases. 
Figure [9] shows that the stochastic hybrid is able to correctly reproduce the rate of exponential 
relaxation of the piston toward equilibrium with many fewer particles than the purely particle runs, 
while the deterministic hybrid fails. We have observed a slight dependence on the exact details of 
the hybrid calculations such as cell size or the width of the particle subdomain; however, in general, 
the stochastic hybrid has shown to be remarkably robust and successful. At the same time, the 
importance of including thermal fluctuations in the continuum subdomain is revealed as for the 
VACF computations in Section [IV C[ 

The relation to the equilibrium VACF computations in Section IV C| is emphasized by computing 
the VACF C(t) = (v x (t)v x (0)) for the piston in its state of true equilibrium 4 , k^T = 1. The 



3 At thermodynamic equilibrium with a common temperature T, the frequency of the thermally-driven oscillations 
can be estimated using a quasi-adiabatic harmonic approximation to be ui 2 ~ NksT/ [Mx(L x — x)\ (see also an 
alternative derivation in Ref. 1691) and the amplitude of the oscillations can be estimated to be on the order of 
Aa: 2 « x(L x — x)/N, where N is the number of fluid particles per chamber. 

4 In Ref. Ref. 69 the autocorrelation function for the piston bead position is used to extract the rate of exponential 



40 



C(t) 

l.OxlO" 3 



7.5x10 



5.0x10" 



2.5x10" 



5x10 



-2.5x10" 



-5.0x10" 



Particle 
Stoch. w =2 



Particle 



Stoch. hybrid 

Det. (w p =4) 
Det. xlO 




50 



100 



150 



200 



250 



Figure 10: The velocity autocorrelation function for a rigid piston of mas M/m = 1000 at thermal equilibrium 
at fceT = 1 (thus, the expected initial value is C(0) = 10~ 3 ), and macro cells of size 4 x 40 micro cells. The 
deterministic hybrid gives much smaller effective temperature T e g of the piston and a negative dip in C (t) 
at short times, but when magnified by an order of magnitude it reveals the correct shape at longer times. 
The inset focuses on the initial decay of the VACF and shows that even increasing the width of the particle 
region to wp — 8 (the whole box has a length of 60 continuum cells) does not help the deterministic hybrid 
much whereas the stochastic hybrid gives the correct initial decay. 



VACF is rather complex due to the interplay of short-time kinetic effects, dissipation, and sound 
reflections from the walls, however, our focus here is to simply compare against the purely particle 
simulations and not to understand all the features in the VACF. The results, shown in Fig. [TO] , 
reveal that the piston does not equilibrate at the correct effective temperature T e g = MC(0)//cb 
in the deterministic hybrid calculations. Notably, just as we found for the massive bead in Section 



IV C the piston has a kinetic energy that is markedly lower (half) than the correct value C(0) = 
ksT/M = when fluctuations are not consistently included in the continuum region. Since 
the quasi-equilibrium temperature of the piston 5 plays a crucial role in all of the kinetic theories 
|68 1 169 1 ITT] f° r the effective heat transfer, it is not surprising that the deterministic hybrid gives the 
wrong answer. What is even more striking is that increasing the width of the particle subdomain wp 



rate toward true equilibrium, however, a direct non-equilibrium calculation as we perform in Fig. |9]is more efficient 
and illustrative for our purpose. 
5 It is predicted that the piston equilibrates at a temperature that is approximately the geometric mean of the left 
and right temperatures [73j . 
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on each side of the piston, as measured in units of continuum cells (Fig. [7]shows the case wp = 2), 
barely improves the accuracy of the VACF at short times, showing that the whole spectrum of 
fluctuations affects the initial rate of decay of the piston. At the same time, the deterministic 
hybrid does give the correct shape of the VACF at longer times, as revealed by magnifying the 
VACF for wp = 4 by an arbitrary factor of 10 to bring it in close agreement with the particle result 
at longer times 6 . This is also not unexpected since the VACF at longer times is dominated by 
mechanical vibrations and dissipation, and is analogous to what we find for the dynamic structure 
factor when it is calculated by a deterministic hybrid scheme. 



V. CONCLUSIONS 



We have described a hybrid particle-continuum algorithm for simulating complex flows and ap- 
plied it to several non-trivial problems. The algorithm couples an ideal stochastic particle fluid al- 
gorithm with a fluctuating hydrodynamic continuum solver using a direct dynamic coupling where 
the continuum solver supplies Dirichlet-like (state) boundary conditions for the particle region, 
while the particle region supplies von Neumann-like (flux) boundary conditions to the continuum 
solver. The continuum solver computes a provisional solution over the whole domain that then 
gets replaced with the particle data in the particle region and also gets corrected (refluxed) in the 
cells bordering the particle subdomain. We described the components necessary to extend previous 
variants of this rather general coupling methodology [20] [2l"1 [29] to use a recently-developed variant 
of the DSMC particle method suitable for dense fluids |18| 1ST] , as well as an explicit conservative 
compressible fluid solver that accurately accounts for thermal fluctuations in the Navier-Stokes 
equations [9] . By turning the fluctuating fluxes off in the continuum solver we can trivially trans- 
form our stochastic hybrid method to a deterministic hybrid method closer to commonly-used 
hybrid schemes. 



In Section IV we applied our stochastic hybrid method to several challenging problems and 
demonstrated that it could obtain the correct answer with significantly fewer particles and thus 
significantly less computational effort than a purely particle simulation. We used purely particle 
runs as a "gold" standard against which we judge the accuracy of the hybrid method, consistent 
with using the hybrid method for situations in which a purely continuum description is unable to 
capture the full physics and a particle method is necessary in some portion of the physical domain. 



6 The magnification factor required to bring the long-time VACF for the deterministic hybrid in agreement with the 
correct result decreases with increasing top, for example, it is ~ 20 for wp — 2 and ~ 4 for wp — 8. 



42 



For example, particle methods such as DSMC are essential to correctly resolve the flow at a high 
Mach number shock [21], or in kinetic flow regions such as the wake behind a fast-moving body 
[20] or small channels in a micro-electromechanical device [S] . 



In Section IV A we showed that there is a small mismatch between the density and temperature 
in the particle and continuum regions, caused by using instantaneous fluctuating values of the 
hydrodynamic variables when generating velocities for the reservoir particles, and that the effect 
is of order Nq 1 , where Nq is the average number of particles in one cell. In Appendix [X] we 
demonstrated that this is not an artifact of the method but rather of an inherent difficulty in 
stochastic methods where instantaneous values are used to estimate means. In section IIVBI we 
calculated the dynamic structure factor using the hybrid method for a quasi two-dimensional 
situation and a large wavevector k that is obliquely incident to the particle-continuum interface, 
and confirmed that spontaneous sound and entropic fluctuations are transmitted correctly through 
the interface. We also calculated the dynamic structure factor for a finite system bounded by two 
adiabatic walls and observed excellent agreement with theoretical calculations given in Appendix 

m 



In Section IV C we studied the diffusive motion of a large spherical neutrally-buoyant bead 
suspended in a fluid of I-DSMC particles in three dimensions by placing a mobile particle subdomain 
around the suspended bead. We computed the velocity autocorrelation function (VACF) for two 
bead sizes using the stochastic and deterministic hybrids as well as purely particle simulations. 
We found that the stochastic hybrid correctly reproduces the VACF over all time scales, while the 
deterministic hybrid under-estimates both the kinetic energy of the bead and the magnitude of the 



tail of the VACF for sufficiently large beads. Finally, in Section IV D we applied the hybrid scheme 
to a non-equilibrium quasi one-dimensional version of the adiabatic piston problem |67[ [68] , a classic 
example of the importance of fluctuations in establishing global thermal equilibrium. The two- 
dimensional particle region was placed around the piston and an event-driven algorithm was used 
to handle the interaction of the fluid with the piston. We again found that the stochastic hybrid 
was able to reproduce the purely particle results correctly, while the deterministic calculations 
under-estimated the relaxation substantially for sufficiently massive rigid pistons. 

Our results for both the VACF and the adiabatic piston clearly demonstrated that a large 
massive suspended body cannot equilibrate at the correct Boltzmann distribution unless thermal 
fluctuations are consistently included in the full domain, including the continuum region in hybrid 
methods, even if a large particle subdomain is used. This points to an increased importance of 
the long-wavelength, and thus slowly- decaying, hydrodynamic fluctuations. Massive and large 
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suspended bodies have longer relaxation times and thus it is not surprising that slower-decaying 
fluctuations play a more prominent role for them than for smaller suspended particles. However, 
a better theoretical understanding of these observations is necessary in order to establish the 
importance of hydrodynamic fluctuations in general. At the same time, our results make it clear 
that fluctuations should be included in the continuum region in hybrid methods, consistent with 
the particle dynamics, rather than treating the fluctuations as "noise" from which the continuum 
solver ought to be shielded. 

Fluctuating hydrodynamics has successfully accounted for thermal fluctuations in a variety of 
problems. At the same time, however, the nonlinear stochastic partial differential Landau-Lifshitz 
Navier-Stokes equations are mathematically ill-defined and require a cutoff length and/or time 
scale to be interpreted in a reasonable sense. The linearized equations can be given a well-defined 
interpretation, however, they are physically unsatisfactory in that they require a base state to lin- 
earize around which may itself be time-dependent and unknown or even be affected by fluctuations, 
as in the adiabatic piston or diffusing shock problems [6] [19]. We have numerically observed that 
using small continuum cells leads to worse results even if the linearized LLNS equations are used, 
thus formally avoiding the difficulties with the increased relative magnitude of the fluctuations. 
This is consistent with the expectation that a continuum description is only applicable on length 
scales and time scales sufficiently larger than the molecular size and molecular collision time. In 
our experience using more than 75 particles per cell leads to a good match between the hybrid and 
particle runs; however, a better theoretical understanding of the proper inclusion of fluctuations in 
hydrodynamics is a necessary future development. 

Our implementation is at present serial and our runs are therefore limited by the CPU-intensive 
collision procedure in the I-DSMC particle algorithm. Some of the examples we presented utilized 
the mixed event- and time-driven particle algorithm developed in Ref. [30]. The event-driven 
component of this algorithm is notoriously difficult to parallelize. However, a purely time-driven 
particle algorithm can easily be parallelized, as can the purely continuum solver. We plan to 
implement a parallel hybrid scheme in the future in order to enable the study of realistic system 
sizes. At the same time, however, reaching long time scales will necessarily require time steps 
beyond the small ones required by particle methods and explicit continuum schemes. 
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Appendix A: FINITE-SIZE MISMATCH BETWEEN PARTICLE AND CONTINUUM 

DESCRIPTIONS 

Consider a simple particle/continuum hybrid consisting of a one-dimensional system with an 
ideal particle fluid on one side and a fluctuating continuum solver on the other. At equilibrium the 
mean density and temperature are po and To, respectively; the mean fluid velocity is taken as zero. 
The problem we want to consider is whether the one-sided fluxes of mass, momentum, and energy 
due to the reservoir particles that enter the continuum subdomain are correct, assuming that the 
means and variances of the hydrodynamic variables in the reservoir cells are correct. 

In the continuum calculation instantaneous density and temperature, p and T, are computed 
and used to generate particles for injection. The one-sided fluxes of mass, momentum, and energy 
are known from kinetic theory to be 



F E (p,T) = J^-pT 3 '*. 
V m d TT 

The correct mean equilibrium mass flux is F p (po, To); however, the mean value of the particle flux 
does not equal this correct value since 

{pT 1 / 2 ) = {p)(T 1 / 2 )^(p)(T) 1 / 2 ; 

note that density and temperature are instantaneously uncorrelated. Similarly, the mean energy 
flux is also incorrect; interestingly the momentum flux is correct since it has a linear dependence 
on temperature. 

To find the errors in the fluxes, we write T = To + ST and for a given power exponent a we have 



(T a ) = (T$)((l + 5T/T ) a )^{TS) 



1 + \a{a - 1) 



(ST 



2\ 



1 



(T$) 



1 + 



a(o-l)( 7 -l) 



2N 

where A^o = poV c /m is the number of particles in a continuum cell (of volume V c ). From this we 
have that, for a monatomic gas (7 = 5/3), 

(F p (p,T)) = Fp (po,T )(l-^ 
{Fj(p,T)) = Fj(p ,T ) 
(F E (p,T)) = F E (p ,T )(l + 



12N 
1 
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This result shows that there is no way for all three mean fluxes to be correct when the particles 
are generated from the Maxwell-Boltzmann distribution using instantaneous, fluctuating values of 
density and temperature. The mismatch is of the order Nq 1 , where Nq is the number of particles 
per macro cell, and therefore the mismatch gets worse as the macroscopic cells become smaller and 
the (relative) fluctuations become larger. As our numerical results show, because of this mismatch 
between the particle and hydrodynamic descriptions it will be impossible for the particle and 
continuum regions to reach a common equilibrium state. 

Appendix B: DYNAMIC STRUCTURE FACTOR WITH ADIABATIC WALLS 

Dynamic structure factors are easily calculated in the bulk by using the spatio-temporal Fourier 
transform of the linearized LLNS equations [53]. For finite domains, such as slab geometries, there 
are results in the literature but they are often restricted to some simplified models or complex 
non-equilibrium situations |74[ [75] . We therefore derive here the equilibrium dynamic structure 
factors for a fluid in-between two adiabatic hard walls by solving the LLNS equations with the 
appropriate boundary conditions. 

Boundary conditions change the Hilbert space in which a solution is to be sought and the cor- 
responding basis functions (eigenfunctions of the generator with the specified BCs). For the LLNS 
equations with adiabatic boundaries (i.e., slip insulating walls) along the x axes, the appropriate 
basis functions are cos{kx) for density and temperature, where k = pir/L is the wavevector enu- 
merated by the wave index p S Z + , and sin(/cx) for the velocities |74| I75j. as compared to the ones 
for "bulk" conditions (periodic boundaries), exp(2qirx/L), with k = 2qn/L and wave index q £ Z 
|53]. For thermal walls (constant-temperature stick boundaries) there does not appear to be a 
simple basis. 

White noise has a trivial expansion in either the sine or cosine basis sets, namely, all of the 
coefficients are i.i.d. Gaussian random variables with mean zero and variance 2. The generator of 
the Navier-Stokes equations separates wavevectors/frequencies into the the same (k,uj) equations 
as for "bulk" (periodic BCs), and therefore the dynamic structure factor, if expressed in the given 
basis set, has the same familiar form \7i\ ITS] . In particular, 



p{x, t) = po + 22 p p (t) cos(pirx/L) = p + / due luJt p P)L u cos(pirx/L) 




where the different p's and u's are uncorrelated 
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where S p>UJ = S k=pn / L ^ denotes the usual bulk dynamic structure factor. 

Now, we need to convert this result to the more usual Fourier basis, exp(2(/7rx/L), since this is 
how dynamic structure factors are defined. From the Fourier inversion formula, and the orthogo- 
nality of the cosine basis functions, we have 

1 f L 1 i °° f L 

p q = — \ dxp(x,t)e~ 2qnx ' L = -p P =2q ~ yVft / dxcos(pnx/L) sm(2qirx/L). 
L Jo 2 L ~{ Jo 

By performing the integration explicitly we get 



giving the dynamic structure factor 



i iiq sr^ 

P q= 2 p p=2q + — 2, 2 _ 42 : 
p odd 



p odd " " 



Each of the terms under the sum gives an additional peak at to = ck = pcn/L, where p is odd, 
which arises physically because of the standing waves that appear due to reflections of the sound 
waves from the walls. Only the first few terms need to be kept to get most of the power (the 
total integral over u is one), and it can be shown (by simple numerical comparison or explicit 
summation) that the above is equivalent to the more opaque Eq. (12) in Ref. [75] (set 7 = for 
no shear). 



Equations identical to (Bl) hold for temperature and the velocity components parallel to the 
wall. For the perpendicular velocity {v x ) a sine basis is more appropriate, and a similar calculation 
gives the dynamic structure factor for the finite system in terms of the bulk one, 

D k=2irq/L,u) o% "i _2 / , 



k=2irq/L,u> 2 k ' u 7T 2 ^ ( r> 2 - 4a 2 ) 2 ' 

p odd ^ " i 



Appendix C: HANDLING OF ADIABATIC AND THERMAL WALLS IN THE 

CONTINUUM SOLVER 



Solving the LLNS equations with non-periodic boundaries requires some special handling of the 
stochastic fluxes at the boundaries, which are assumed to coincide with faces of the continuum 
grid. As discussed in Refs. [HI |76], the numerical discretization of the Laplacian operator L, the 
divergence operator D, and the gradient operator G should satisfy a discrete fluctuation-dissipation 
balance condition L = DCG = —DCD*, where C is a dimensionless covariance matrix for the 
stochastic fluxes that are generated using a random number generator on each face of the grid. For 
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one-dimension with periodic boundaries it is well known that the standard face-to-cell divergence, 
cell-to-face gradient, and three-point Laplacian second-order stencils satisfy L = DG and thus 
C = I (the identity matrix) works and it is in fact trivial to implement algorithmically ^\ 176] . 
When boundaries are present, the stencils near the boundaries are modified to take into account 
the boundary conditions. 

Algorithmically, ghost cells extending beyond the boundaries are used to implement modified 
finite-difference stencils near the boundaries. The numerical scheme continues to use standard 
divergence D (face-to-cell) and gradient (cell-to-face) G = —D* stencils but implements a modified 
Laplacian operator due to special handling of the ghost cells. If a Dirichlet condition is imposed on 
a given variable (e.g., a fixed wall temperature), then the ghost cell value is obtained by a linear 
extrapolation of the value in the neighboring interior cell (inverse reflection). If a von Neumann 
condition is imposed, on the other hand, then the ghost cell value is set equal to the value in the 
neighboring interior cell (reflection). This gives discrete operators that can be represented by the 
following banded matrices near the left corner (first cell) of a one-dimensional domain 



D 



-1 1 ••• 
-1 1 ... 
0-1 ••• 



-(2 -a) 1 ••• 
1 -2 1... 
1 -2 ••• 



where a = — 1 for a Dirichlet condition and a = 1 for a von Neumann condition. It is easy to verify 
that L = —DCD* is satisfied with the following diagonal scaling matrix 



C 



(3 
1 
1 



where (3 = 1 — a. 

This direct computation shows that in order to satisfy the discrete fluctuation-dissipation bal- 
ance condition the diagonal element of C corresponding to the cell face that touches the boundary 
ought to be set to 2 for a Dirichlet and to for a von Neumann condition. This means that 
the corresponding component of the stochastic flux needs to be generated using a random normal 
variate of variance 2 for Dirichlet, and set to zero for a von Neumann condition. 

Finally, for density, the ghost cell values are generated so that the pressure in the ghost cells 
is equal to the pressure in the neighboring interior cell, which ensures that there is no unphysical 
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pressure gradient in the momentum equation across the interface. There is also no stochastic mass 
flux through faces on the boundary independent of the type of boundary condition at the wall. 
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